.. _tiskitpy.ComplianceNoise_example:

==============================
SeafloorSynthetic example code
==============================

.. code-block:: python

    """
    Create synthetic OBS training data

    - Read in real data from a quiet continental site
    - Create a noise model using the SeafloorSynthetic class with default values
    - Add the two together
    - Save the result
    - Also save:
        - The original input data
        - The calculated compliance
    
    """
    from pathlib import Path

    import numpy as np
    from obspy import read, read_inventory, UTCDateTime
    import matplotlib.pyplot as plt

    from tiskitpy import SpectralDensity, SeafloorSynthetic, PSDVals

    data_file = 'data/G.TAM_2010059-2010069.mseed'  # post-Maule eq
    inv_file = 'data/G.TAM.2010.station.xml'
    wdepth = 2400
    station = 'SYNV1'
    wt = 'prol4pi'  # hamming and prol1pi don't have enough broadband noise rejection

    # Read the real data and its metadata
    real_data = read(data_file, 'MSEED')
    resp_trace = real_data[0].copy()

    # Change the start time to "hide" its provenance
    resp_trace.stats.start_time = UTCDateTime(2024,1,1)
    inv = read_inventory(inv_file, 'STATIONXML')
    # Extract the instrument response
    s_response=inv.select(channel='LHZ')[0][0][0].response

    # Create the noise model and synthetic data stream
    noise_tilt_max = PSDVals.sloped_freqs_and_values(-180, -30, -3, 0.1, .25)
    noise_model = SeafloorSynthetic(noise_tilt_max=noise_tilt_max,
                                  noise_tilt_variance=30)

    # PLOT noise model using SeafloorSynthetic's intrinsic method
    noise_model.plot(outfile='noise_model.png')

.. image:: images/12_PeridocTransient_timing_clips_1000.png
   :width: 564
   
.. code-block:: python

    # Create the synthetic and sources data streams
    data_synth, sources, inv_synth = noise_model.streams(
        real_data[0], s_response=s_response, p_response = 100.,
        station=station)
    noisePSDs = noise_model.PSDs

    # Compare the PSDs of the noise sources' waveforms to their theoretical valus
    sd_sources = SpectralDensity.from_stream(sources, windowtype=wt)
    fig, ax = plt.subplots()
    for component, color in zip(('IGZ', 'NOS', 'NTZ_max', 'NTZ_min'), ('r', 'b', 'g', 'g')):
        ax.semilogx(noisePSDs[component].freqs, noisePSDs[component].values, color, label=component)
        ax.semilogx(sd_sources.freqs, 10*np.log10(sd_sources.autospect('*' + component[:3])), color+'--')
    ax.set_ylim(-200, -80)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('dB ref 1 m/s^2/sqrt(Hz)')
    plt.suptitle(f'Z theoretical versus sd(stream, {wt}) noise levels')
    plt.legend()
    plt.show()

.. image:: images/11_ComplianceNoise_Zsources.png
   :width: 564
   
.. code-block:: python

    # PLOT synthetic time series
    data_synth.plot(equal_scale=False)


.. image:: images/11_ComplianceNoise_Synth_Stream.png
   :width: 564
   
.. code-block:: python

    # PLOT synthetic PSD versus SeafloorSynthetic components
    sd_synth = SpectralDensity.from_stream(data_synth, inv=inv_synth, windowtype=wt)
    ax = sd_synth.plot_one_autospectra(f'XX.{station}.00.LHZ')
    ylim = ax.get_ylim()
    for label, color in zip(('IGZ', 'NTZ_min', 'NTZ_max', 'NOS'), ('r', 'm', 'c', 'g')):
        ax.plot(noisePSDs[label].freqs, noisePSDs[label].values,
                color=color, label=label, linestyle='-.')
    ax.set_ylim(ylim)
    ax.legend()
    plt.gcf().suptitle('Synthetic Data PSD and its sources')
    plt.show()

.. image:: images/11_ComplianceNoise_SynthZ_PSDs.png
   :width: 564
   

.. code-block:: python

    # Add the real and synthetic data together
    data = data_synth.copy()
    data.select(channel='LHZ')[0].data += real_data.select(channel='LHZ')[0].data
    data.select(channel='LH1')[0].data += real_data.select(channel='LHN')[0].data
    data.select(channel='LH2')[0].data += real_data.select(channel='LHE')[0].data

    # PLOT synthetic + real time series
    data.plot(equal_scale=False)


.. image:: images/11_ComplianceNoise_SynthReal_Stream.png
   :width: 564
   
.. code-block:: python

    # PLOT synthetic + real PSD versus SeafloorSynthetic components
    sd_synth = SpectralDensity.from_stream(data, inv=inv_synth, windowtype=wt)
    ax = sd_synth.plot_one_autospectra(f'XX.{station}.00.LHZ')
    ylim = ax.get_ylim()
    for label, color in zip(('IGZ', 'NTZ_min', 'NTZ_max', 'NOS'), ('r', 'm', 'c', 'g')):
        ax.plot(noisePSDs[label].freqs, noisePSDs[label].values,
                color=color, label=label, linestyle='-.')
    ax.set_ylim(ylim)
    ax.legend()
    plt.gcf().suptitle('Synthetic + Real')
    plt.show()

.. image:: images/11_ComplianceNoise_SynthReal_PSDs.png
   :width: 564

.. code-block:: python

    sd_data = SpectralDensity.from_stream(data, windowtype=wt)
    # PLOT synthetic+real PSD
    sd_data.plot()

.. image:: images/11_ComplianceNoise_SynthReal_PSD.png
   :width: 564

.. code-block:: python

    # PLOT synthetic+real coherence
    sd_data.plot_coherences()

.. image:: images/11_ComplianceNoise_SynthReal_Coherences.png
   :width: 564